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ABSTRACT 



We apply linear response theory to a general, inhomogeneous, stationary stellar 
system, with particular emphasis on dissipative processes analogous to Landau 
damping. Assuming only that the response is causal, we show that the irreversible 
work done by an external perturber is described by the anti-Hermitian part of 
a linear response operator, and damping of collective modes is described by the 
anti-Hermitian part of a related polarization operator. We derive an exact formal 
expression for the response operator, which is the classical analog of a well-known 
result in quantum statistical physics. When the self-gravity of the response can 
be ignored, and the ensemble-averaged gravitational potential is integrable, the 
expressions for the mode energy, damping rate, and polarization operator reduce 
to well-known formulae derived from perturbation theory in action-angle vari- 
ables. In this approximation, dissipation occurs only via resonant interaction 
with stellar orbits or collective modes. For stellar systems in thermal equilib- 
rium, the anti-Hermitian part of the response operator is directly related to 
the correlation function of the fluctuations. Thus dissipative properties of the 
system are completely determined by the spectrum of density fluctuations — the 
fluctuation-dissipation theorem. In particular, we express the coefficient of dy- 
namical friction for an orbiting test particle in terms of the fluctuation spectrum; 
this reduces to the known Chandrasekhar formula in the restrictive case of an 
infinite homogeneous system with a Maxwellian velocity distribution. 
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1 Introduction 

The aim of this paper is to apply to stellar dynamics two of the powerful tools that have been 
developed in other branches of statistical physics: linear response theory and the fluctuation- 
dissipation theorem. 

The apphcation of the methods of statistical physics to self-gravitating stellar systems is 
challenging for several reasons: the interparticle forces are long-range, the systems are intrin- 
sically inhomogeneous, and true thermodynamic equilibrium does not exist (e.g. Lynden-Bell 
& Wood 1968). Thus, in practice, most of our understanding of relaxation in stellar sys- 
tems comes from analyses based on the following related approximations (Jeans 1913, 1916; 
Chandrasekhar 1942; Binney & Tremaine 1987; Spitzer 1987): 

1. Local approximation: The stellar system is assumed to be infinite and homogeneous, 
and the force field from the equilibrium system is neglected; thus the unperturbed 
stellar orbits are straight lines at constant velocity. The validity of this approximation 
requires that most of the relaxation arises from close encounters, which have impact 
parameter b <^ R where R is the size of the system. The failure of the local approxi- 
mation at large impact parameter is why the Coulomb logarithm is poorly determined 
in stellar systems. 

2. Markov approximation: The interactions between stars are treated as sequential binary 
encounters of negligible duration, so that the evolution of the one-particle distribution 
function (hereafter df) can be treated as a Markov process. This is a reasonable ap- 
proximation for short-range forces such as those in neutral gases, but more questionable 
in a stellar system because of the long range of the gravitational force. The Markov ap- 
proximation is closely related to the local approximation because only local encounters 
can plausibly be considered to have negligible duration. 

3. Diffusion approximation: Most of the relaxation is assumed to arise from weak encoun- 
ters, in which the orbital deflection is small; in practice this requires b ^ Bq = Gm/a'^ 
where m and a are a typical stellar mass and velocity. In this approximation, the evolu- 
tion of the stellar orbits is a diffusion process that can be described by a Fokker-Planck 
equation. The diffusion approximation is consistent with the local approximation if 
bo -C R, so that most encounters are weak (6 ^ bo) while still local (6 -C i?); in this 
case the contribution to the relaxation from encounters with b < bo is smaller than the 
contribution from & > 6o by of order the Coulomb logarithm In A = In it!/6o- With the 
approximations listed so far, the effects of relaxation on a stellar orbit are completely 
described by a dynamical friction force, which gives the mean change in velocity per 
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unit time (Afj), and a diffusion tensor, which gives the mean-square change in velocity 
per unit time {AviAvj). 

4. Neglect of self-gravity: The self-gravity of the response of the stellar system is ne- 
glected. Thus, for example, dynamical friction arises from the gravitational drag 
exerted on a body by the wake it creates, but most calculations of dynamical fric- 
tion neglect the effects of the self-gravity of the wake in determining its amplitude and 
shape. The error introduced by this approximation can be significant: Weinberg (1989) 
and Hernquist & Weinberg (1989) find that including the self-gravity of the wake can 
suppress dynamical friction on an orbiting satellite by a factor of 2-3. Errors of this 
magnitude are often excused because the precise value of the Coulomb logarithm is not 
known in any case; however, if the orbital frequency of the perturbing body is close 
to the frequency of a collective mode of the system, then the effects of self-gravity are 
increased dramatically (e.g. Weinberg 1993). 

The goal of this paper is to summarize some of the insights into fluctuations and dissi- 
pation in stellar systems that can be derived without making any of these approximations. 
Our tools will be linear response theory and the fluctuation-dissipation theorem (e.g. Martin 
1968, Forster 1975, Landau & Lifshitz 1980, Reichl 1980, Sitenko 1982, Klimontovich 1986, 
Kubo et al. 1991). 

The intimate relation between fluctuation and dissipation in stellar systems was first 
recognized by Chandrasekhar (1943), who argued that a star's random walk in velocity due 
to stochastic gravitational forces must be balanced by a drag force (dynamical friction) if the 
stochastic process is to leave a Maxwellian distribution invariant. The most general form 
of this relation is described by the fluctuation-dissipation theorem. This theorem relates 
the fluctuations in a dynamical system in thermal equilibrium, described by the correlation 
function, to the rate at which the system absorbs energy from a weak external field, described 
by the response operator. In the context of stellar systems, the theorem relates the dynamical 
friction force to the diffusion tensor (c/. eq. |151| below); however, we shall see that this 
relation is far more general than the usual derivation of either the dynamical friction force 
or the diffusion tensor individually. 

The conceptual importance of the relation between dynamical friction and stochastic 
forces was recognized by Bekenstein & Maoz (1992) and Maoz (1993), who concentrated 
on providing a unified derivation of both effects in a homogeneous system; in contrast, we 
shall focus on the general properties of drag and fluctuations that can be derived without 
reference to specific systems. 

We begin §0 with a review of linear response theory applied to a stationary inhomoge- 



-4- 



neous system, stressing the constraints that causahty places on the analytic properties of the 
response operator; for example, we derive the Kramers-Kronig relations between the Her- 
mitian and anti-Hermitian parts of this operator and show that dissipation in the system is 
determined by the anti-Hermitian part of the response operator. We also define the polariza- 
tion and dielectric operators, and derive their relation to the response operator. We employ 
these operators to examine collective modes of the system in deriving an expression for 
the energy of a mode in terms of the Hermitian part of the polarization operator, and show 
that its (Landau) damping rate is determined by the anti-Hermitian part of the polarization 
operator. In §^ we derive an exact formal expression for the response operator in a general 
Hamiltonian system, which has a direct analog in quantum systems. For systems described 
by integrable potentials, we give an explicit form for the polarization operator in terms of 
the one-particle df. We describe fluctuations in the system using the correlation function in 
deriving the symmetries that are imposed on the correlation function by the principle of 
microscopic reversibility. In §^ we prove the fluctuation-dissipation theorem, which relates 
the response operator and the correlation function for an isothermal system. We describe 
fluctuations in non-isothermal systems and the dressed-particle approximation in §^ In §H 
we apply these results to derive general relations between dynamical friction and relaxation 
in isothermal stellar systems. We end with a discussion in §P, and a summary of the main 



formulae and results in 510 



Throughout this paper, we shall apply the term "equilibrium" to a stellar system when 
the one-particle df is a solution of the time-independent collisionless Boltzmann equation; 
such a system is stationary except for slow evolution due to relaxation. We reserve the term 
"thermal equilibrium" for a system whose A^-particle df is an exponential function of the A^- 
particle Hamiltonian. Most of our results apply to any equilibrium stellar system, although 
some apply only to time-reversible systems (in which the Hamiltonian and the equilibrium 
one-particle df are invariant under time- reversal), to systems in which the Hamiltonian is 
integrable, or to systems in thermal equilibrium. 



2 Linear response theory 



In this section we examine the response of an equilibrium stellar system to small potential 
perturbations, assuming that the induced density depends linearly on the strength of the 
perturbation through an — as yet unspecified — response operator. It is often difficult or 
impossible to find an analytic expression for the response operator. Nevertheless, many of 
its properties follow directly from causality and linearity. 
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2.1 The response operator 



We consider an equilibrium stellar system that is subjected to a small perturbation $e(r, t) 
from an external potential. The density perturbation induced in the system, ps{r,t), can be 
expressed in terms of a linear response operator R{r, r', r), defined by 



Ps{r,t) = J dr'dt'R{r,r',t-t')<l>e{r',t'y, 



causality requires R{r, r', r) = for r < 0. We stress that does not include the gravita- 
tional potential arising from the response density ps] in this respect the response operator 
can be contrasted with the polarization operator defined in §|2.2| below. 



In this section we consider the general properties of the response operator, deferring the 
derivation of explicit forms for R{r, r', r) until §^ We start by taking Fourier transforms, 
which we denote by replacing the variable r by u, e.g. $e(i'5 t) = duj^e{T', u) exp(— itut). 
Equation (|1|) simplifies to 

Ps{t, uj) = 27r J dr'Riv, r', u;)$e(r , u). (2) 

Since -R(r, r', r) is real, we must have 

R*{t,v\uj) = R{y,v\-uj). (3) 

We can analytically continue the response operator so that it is defined for complex frequen- 
cies z {cf. eq. |T8| with u z). Causality requires that R{r, r', z) has no poles in the upper 
half plane. 

Equation (||) can be written as an operator equation, 

Ps = 27rR(c^)$e- (4) 

If we define the inner product as 

{i:,(P)= j dvrij)(t){v)dv, (5) 
then the adjoint operator Yi^uj) satisfies 

(RV, 0) = (V^, R0) for all V(r), 0(r), (6) 

and is given by 

R\v,v',u:) = R*{v',v,u). (7) 
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We can always write the response operator in the form 

R{r, r', u) = Rnir, r', u) + i?A(r, r', u), (8) 

defined by 

RHir,r',uj) = ^Rir,r\uj) + lR*{r',r,u), 

Ra{v,v\uj) = ii?(r,r',a;)-|/?*(r',r,cu). (9) 

The operators Rh and R^ satisfy the relations 

R^H{r,r',cu) = RUr',r,u) = RHir,r',uj), R\{v,t' .u) = R\{v\t,uj) = -Ra{v,v\u:). 

(10) 

Thus Rh is Hermitian and Ra is anti-Hermitian. Using equations (|^) and @ it is easy to 
show that 

R*H{r,Y\uj) = Rh{y,v' ,-uj), i?^(r, r',u;) = i?A(r,r', -cj). (11) 

In terms of time rather than frequency, we may write 

i?(r, r', r) = R,{v, r', r) + /?„(r, r', r), (12) 

where 

i?e(r,r',r) = J dujRH{r,r',uj)e-'''\ i?o(r,r',r) = j duRAij^r' ,u)e-'''\ (13) 

It is straightforward to show using (pJ]) that Re and Ro are real. Using ([TOD we can show 
that 

/2e(r,r',r) = i?e(r,r,-r), i?„(r, r', r) = -i?„(r', r, -r). (14) 
Moreover, since i?(r, r', — r) = for r > we must have 

i?e(r, r , r) = i?„(r, r', r), i?e(r, r', -r) = -R^iv, r', -r), r > 0. (15) 

If -R(r, r', r) is well-behaved then -R(r, r', z) ^ as Im(2;) — > cxd, where z is the complex 
frequency. Since -R(r, r', z) has no poles in the upper half plane, we can write 

du' ^ ' ' ' = 0, 16 

3 u' ~ UJ + Zl] 

where t] > 0, since we can close the contour at Im(2;) = +oo. We next make use of the 
identity 

li"^ -—— = 'P {——] -'i''^^S'^iv)Six -y), (17) 

v^o X — y + IT] \x — y J 
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where V denotes the Cauchy principal value and 5 denotes the Dirac delta function. Then 
we have 

R{r, r', u;) = / du;'R{r, r', u')V , (18) 

Equating the Hermitian and anti-Hermitian components, 

n J-oo \uj' — u J 

RA{r,r',uj) = r duj'RH{v,v\uj')v(^^). (19) 

These are the Kramers-Kronig relations, which follow from causahty and require no other 
assumptions about the dynamics of the stellar system. They can be analytically continued 
to complex frequencies z — uj -\- ir]. 

Another relation is obtained by examining the integral 

i?(r, r', z)dz, (20) 



z"^ + 



where s is real and positive, and the integral is taken along the real axis and closed in the 
upper half-plane (e.g.. Landau & Lifshitz 1980). Since i?(r, r', 2;) has no poles in the upper 
half-plane, the only contribution to the integral comes from the pole at 2; = is; thus 

fl(r.r-...) = -ir -''^^'/f )^ . (21) 

TT J-oo o;^ -|- 

Equating the Hermitian and anti-Hermitian components, 

RH{r,r',is) — — — / da;^"'^'^^" ' " ' 

TT J-: 



, uRa^y, r', u) 
up' -|- 



-00 



RA{r,r,zs) = / duj r— ^ — . 22 

TT J-oo -l- s 



2.2 The polarization and dielectric operators 



It is sometimes useful to introduce a different measure of the linear response: the polar- 
ization operator relates the induced density to the total potential produced by an external 
perturbation, — + 



Ps{r,t)^ j dv'dt'P{v,v',t-t')^t{v',t'), 



(23) 
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where $s is associated with the induced density ps by Poisson's equation, V^$s = inGps- In 
frequency space P(r, r', u) has the same analytic properties described in § pTl| for _R(r, r', u); 
for example, P(r, r', u) satisfies Kramers-Kronig relations. In general P is easier to compute 
explicitly than R because it does not depend on the self-gravity of the response, which is 
already included in $j. 

It is straightforward to show that the two operators are related by the following nonlinear 
integral equations, 

R{r, r', u) = P(r, r', to) + 27r J dxdx'P{r, x, cj)^(x, x')i?(x', r', u), 

= P(r,r',cj) + 27ry"rfxrfx'i?(r,x,cu)^(x,x')P(x',r',w), (24) 

where 

^(x,x') = -^^ (25) 

is the Coulomb interaction potential. These can be written in operator notation as 

R{z) = P{z) + 27rP(2)*R(^) = P{z) + 27TK{z)^-p (z) . (26) 

Note that * is self-adjoint, = 

Poisson's equation may be written 

<^s{r,uj) = J c/x^(r,x)p,(x,cj) = '^tt J c/xcir'^(r, x)P(x, r', u;)$i(r', cj). (27) 

Thus the external potential is related to the total potential by 

$e(r, cu) = J dv'D{v, r', a;)$i(r', u), (28) 

where 

D{r, r', uj) = 6{r - r') - 2n J dx^(r, x)P(x, r', u) (29) 

is the dielectric response operator. The name arises by analogy to the dielectric constant, 
which similarly relates the total electric field inside a dielectric to the field generated by 
external sources. In operator notation 

$e = D(z)$t where D(z) = I - 27r*P(^). (30) 

The inverse operator satisfies 

^t = 'D-\z)<^e, D-^(^) = I + 27r*R(2). (31) 
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This is easily verified by computing DD ^ and using equation (]26|). The density that gen- 
erates the external potential is related to the total density by 

Pt = A{z)pe where A{z) = I + 27rR(z)*, (32) 

and $e = ^Pe- Its inverse is 

Pe = A~\z)pt where A-\z) = I - 27iP{z)^ . (33) 



The operators D ^{z) and A{z) are singular when the response operator Ii.{z) is singular, 
which occurs at the frequencies of the collective modes of the stellar system (§^. 

Two further identities can be derived using (p6D, 

P(z) = R{z)-D{z), (34) 
K{z) = A{z)-P{z). (35) 



The operators R, P, D and A generally do not commute, and are related through 
the nonlinear integral equations above. These expressions simplify, however, for an infinite 
homogeneous stellar system. In this case, the operators are translationally invariant: the 
spatial dependence of D[r, r', uj), for example, enters only through r — r'. Thus we can define 
a spatial Fourier transform, 

D(k, uj) = J d{v - v')D{y - r', uj)e-'^-^''-'''\ (36) 

so that operating with R, P, D or A reduces to scalar multiplication by a function of the 
wavelength k and the wavenumber uj. In particular for the gravitational Coulomb potential 
given by (|2|), 

*(k) = -'-f. (37) 

Then, 

D(k,c^) = l + ^P(k,w), (38) 

Similar expressions for an infinite homogeneous electron plasma follow with the substitution 
Gm? —y — e^, where — e is the charge of the electron. 
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2.3 Work done by an external potential 



We gain insight into the response operator by considering the work done on the stellar system 
by an external potential $e(r,t); we assume for simplicity that $e ^ as t ^ ±oo. The 
rate at which the external potential does work on a unit mass is — V$e ■ v, where v is the 
velocity; thus the rate of doing work on the stellar system is 

E = - j drV^e ■ pv, 
= y"rfr$eV-(pv), 

= - I dr^^ir,tf-E^, (41) 
where the last line follows from the continuity equation. Thus the total energy change is 

AE, = lEdt = -ldvl dt$,(r,t)^^M = |rfr| rftp(r,t)^5^. (42) 
Using the relation 

/oo roo 
A*{t)B{t)dt = 2n / A*{Lj)B{uj)du (43) 

and the definition of the response operator (|^) we have 

AEs = 27ri J drujdup{r,uj)^l{r,uj) 

= {27ryi J drdr'cudiuR{r,r',uj)<!)^{r',uj)^l{r,uj). (44) 

The contribution of the Hermitian component of the response operator to this integral van- 
ishes, so that 

drdr' / uduRAij, r', cj)<l>e(r', cj)<l>*(r, u) 

J — OO 

//•OO 
drdr' J ^ u;rfu;i?A(r, r', u;)$e(r', c<;)$*(r, cu). (45) 

An important special case of this formula occurs when the external potential is nearly 
monochromatic, with frequency ujq. Then we may write 



/roo 
drdr'RA{r,r,uJo) J duj^e{r' ,uj)^l{r,uj) 

/roo 
drdr'RA{r,r',uJo) / c/t$e(r', t)<l>*(r, t). (46) 
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For example, suppose that the external potential has the form 

$e(r,t) = (7(t)Re[0e(r)e-''^°*], (47) 

where the amplitude g{t) is assumed to be real, vanishingly small in the distant past and 
future, and varies slowly in the sense that \g/g\ = 0(e) <^ \ujo\. Then the component of $e 
with frequency near +Uq is |(7(t)0e(r) exp{—iuQt) , and thus 

AEs = mooo J dtg\t){<f),, R^^e)- (48) 

It is worthwhile to re-derive this result another way. We begin with the last line of 
equation (|4l|). The induced density may be written p(r,t) = Re[ps(r, t)e^'"^"*], where 

/roc 
dr'J^ dTR{r,r',T)Mr')g{t-T)e'^°\ (49) 

Since g{t) changes slowly, we can expand g(t — t) = g(t) — rgit) + O(e^) in (^91) , so that 

p,(r, t) = 2rtg{t) J dv'R{v, r', uJo)<pe{r') + 2mg{t) J dv'^{v, r', cuo)0e(r') + O^e"). (50) 

Thus dp/dt = Re[{-iuJopi + p2)e-''^°^] + 0{e^), where 



dr'R{r, r', ujo)Mr'), P2(r, t) = 2'Kg{t) J dr'-^{r, r', ^o)0e(r'). 

(51) 

We average the rate of doing work over one cycle of the external potential, denoting this 
average by (■). The time average of the integrand in ( ^T] ) can be written 

(^^e^'j = -^(0:Pi - 0ePn + 1{€P2 + 0eP;). (52) 

Finally, writing R in terms of its Hermitian and anti-Hermitian components, using equations 
(|ig) and (|T]), we find 

{E)=W + E,^u (53) 

where 

W = niujQ j (ir(irV*(r)i?^(r,r',cJo)0e(r') = 7rzcJo(0e, Ra^e), (54) 

' ] , (55) 



Eint = -\j drdr'(f)l (r) ^^^^ (r, uJo)(f)e{r') = -| (^(, 
and g{t) is taken to be unity. 



dto 



The quantity W is related to the energy change derived earlier through AE^ = J Wdt; W 
involves the anti-Hermitian part of the response operator and is present even if the external 
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potential is maintained at constant amplitude {g{t) =constant). It represents the rate of 
doing work on the stellar system that is required to maintain the periodic potential; thus 
the anti-Hermitian response Ra is associated with energy absorption (dissipation) by the 
stellar systemQ. In collisionless stellar systems with integrable potentials, the work done by 
the external potential is absorbed by particle resonances or collective modes with frequencies 
near ljq {cf. § [4 .31) . 

The term i^int in ( p3D represents the work required to build up the response to the 
time-varying field. This work must be done even when the dissipation associated with Ra is 
small or zero. We can thus regard E^^t as the total energy associated with the interaction of 
the external field and the stellar system. If the response operator Rh is non-singular at the 
perturbing frequency, this energy will be recovered by the external perturber if the potential 
is turned off adiabatically. 



3 Collective modes 

In the previous section we have considered the induced density response of a stellar system 
that is subjected to an external potential. However, a stellar system may also support a 
self- induced response even when no external perturbation is present {ps ^ even when 
$e = 0). Such a response is called a collective mode; equation (^) implies that the potential 
associated with a collective mode satisfies the dispersion relation 

J dr'D{r, r', z)$,(r', z) = or D(z)$, = 0, (56) 

where z is the complex frequency. The response operator R is singular at the eigenfrequency 
of a collective mode (c/. eq. | and Kalnajs 1971). 

For an infinite, homogeneous medium the collective modes are plane waves, $s oc exp(zk- 
r), and equation (|56|) reduces to an algebraic dispersion relation, D(k,z) = 0. 

3.1 The mode energy 



"'^We refer to this process as energy absorption or dissipation; however, in some coUisionless systems energy 
can be emitted rather than absorbed by this process. 
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We can use the results of § p.3| to determine the energy of a collective mode in the case where 
the imaginary component of the eigenfrequency is small. We imagine perturbing the stellar 
system with an external potential of the form (|17|), where Uq is close to the real part of the 
eigenfrequency of the mode. Averaging equation (5T) over one period 2ti/ujo and writing the 
external potential in terms of the total and response potentials, = — we find 

<-)-/-M>-/-{*-l>-/-{*'l>- 

The first term can be written dUm/dt where 

U^ = ^Jdr {p^,) = ^ [(p„ 0,) + (0„ p,)] = ^ [(P0i, 0,) + (0„ P0i)] (58) 

is the potential energy of the induced disturbance. The second term can be evaluated along 
the lines of equations (P5D-(|55D: 

- Jdr {^t^^ = Wm + Km, (59) 

where 

W.m = TTiuo J drdr'(f)t{r)PAir,r',uJo)(j)t{r') = niuJo{(f)t,PA(l)t), (60) 

Km = -lJ drdr'cp;ir)^iry,uo)Mr') = {<Pu ^0*) • (61) 

Since ujq is close to the eigenfrequency of the mode, we expect that the response of the system 
is strong, ^ \^e\- Thus to a good approximation we can replace ipt by (ps = (pt — 4>e, 
and we find that the total energy associated with a collective mode is 

Em = Km + Um = ^0s, ~^f~^^ + f (0*' ^^(ps) = "^^0 (^4>s, ' (^2) 

and the energy is dissipated at a rate 

Wm = 7iiuJo{(ps, Pa(Ps)- (63) 

Since Um is the potential energy and Em is the total energy, we may identify Km as the 
kinetic energy of the mode. 

The quantity Wm represents the rate of absorption of energy from the mode by the 
stellar system. The amplitude of the mode varies as exp(?7t), where 77 is the imaginary part 
of the eigenfrequency; thus 

_ 1 dEm _ Wm _ i{4>s,PA4>s) 

~ 2Em dt ~ 2Em ~ (0„ Ph,M ' ^ ' 
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where Ph,ui = dPn/duj. These results are only valid in the limit of weak damping, |?7| ^ 
since we have assumed that the eigenfrequency of the mode is close to the real frequency ujq. 

Equation (^) can also be derived directly from the dispersion relation. Taking the inner 
product of ( ^61) with the density of the mode gives 

(p„D(zo)0.) = 0, (65) 

where zq = ujq + irj is the complex eigenfrequency of the mode. We write D in terms of its 
Hermitian and anti-Hermitian parts, and assume that r] is small so that we can expand to 
first order in r], 

(ps, T)H{uJo)(f)s) + {ps, D^(t^o)0s) + ivips, Dh,U^o)'Ps) = 0. (66) 

It is straightforward to show that the first term is real and the remaining terms are purely 
imaginary. Equating the imaginary parts to zero and solving for rj, we recover (|6^) . 



4 The response operator for a Hamiltonian system 



We now examine the dynamics of the stellar system in more detail. We consider a system 
composed of stars, with phase-space coordinates Zj = (rj,Vj), i = 1,...,N (note that 
we define phase space using velocity, not momentum). For simplicity we assume that all 
of the stars are identical, with mass m, although the results we derive can be generalized 
to a range of stellar masses. We denote the coordinates and Hamiltonian of the system in 
6A^-dimensional phase space by Z = (zi, . . . , z^r) and mifo(Z), where 

N 

HoiZ) = 1 ^ + vl/(r, - r,); (67) 

here \l/ is given by equation (|25|) . 

We consider an ensemble of stellar systems, described by a A^-particle distribution func- 
tion (df) /(Z, t), where / is a symmetric function of the variables Zi, . . . , z^y, normalized 
so that / f{Z,t)dZ = 1. Thus the ensemble average of any phase function u{Z) may be 
written (u) = J u{Z)f{Z,t)dZ. The evolution of the df is described by Liouville's equation 

^ + [f,Ho] = 0, (68) 



where [ , ] is a Poisson bracket. 
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The trajectory of the system Z(t) satisfies Hamihon's equations 

f = [Z,i^o]z- (69) 

We shall use Zt- as shorthand for the image of Z after time r under the Hamiltonian fiow 
(|69|); thus, if Z(t) is a trajectory and Z = Z(to), then Z^ = Z(to + t). 

We may derive a generalized response operator for two arbitrary phase functions Us(Z), 
z;s(Z), where s is a parameter that labels these functions. Imagine that we apply a small 
external perturbation mHi{Z, t) to the stellar system. We suppose that mHi can be written 
in the form 

mHi{Z,t) = J ds'v^>{Z)X{s',t). (70) 
The induced perturbation in the expectation of the phase function Us may then be written 

Mi,(t) = j ds'dt'Ruy{s, s', t - t')X{s', t'), (71) 

where i?„t,(s,s',t — f) is a generalized response operator (c/. eq. |1]). 

The perturbation to the df induced by Hi is written /i(Z, t) and satisfies the linearized 
Liouville equation, 

where d/dt denotes the Lagrangian derivative along the unperturbed trajectory. The formal 
solution to this equation is 

/i(Z,t) = -/ rfr[/o,ifi(t-r)]z_,. (73) 
Jo 

Thus we may write 

uis(t) = J dZfi{Z,t)u,{Z) 

//•OO 
dZUsiZ) dT[fo,Hi{t-T)]z_^ 

= -- I dZu^iZ) f ds' rdTX{s',t-T)[fo,v,>]z_^. (74) 
m J J Jo 

Thus the response operator is 

Ruv{s,s',t) = / dZu,{Z)[fo,v,,]z.^, (75) 

m J 
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where 9 is the step function. Using the identity A[_B,C] = [Ai?, C] — -B[yl, C] this can be 
simphfied to 

= ^ f dZfoiZ)[u,iZ,),vAZ)] 

771/ J 

= ^{[u,{t + T),v,,{t)]). (76) 
That is, the response operator is just the ensemble-averaged Poisson bracket of Usit + r) and 

Vs'{t). 

The density-density response operator is defined by setting X(r,t) = $e(r, where 
is the external perturbing potential, and Urit) = fr(t) = p(r, t) where 

N 

p(r,t)=m^5[r,(t)-r], (77) 

i=l 

is the exact density distribution for the A^-body system with coordinates Zj(t) = (rj(t), ri(t)). 
Thus 

i?(r, r', r) = ^([p(r, t + r), p(r', t)]). (78) 

The Poisson bracket is taken with respect to the phase-space coordinates Z = (zi, . . . , zat) 
at time t, which determine ri{t) and rj(t + r) and hence imphcitly determine p(r, t + r) and 
p(r',t) through (0). 

A nearly identical result to (|76D was first derived by Kubo (1957) in the context of 
quantum statistical mechanics. In this case the A^-body df is replaced by the quantum 
mechanical density operator, which satisfies an equation of motion very similar to Liouville's 
equation (|68D ; the primary difference being that the Poisson bracket is replaced by a quantum 
mechanical commutator. 

The expression (|7^) for the response operator is formally exact, but difficult or impos- 
sible to evaluate in practice for realistic stellar systems. However, the analogous expression 
for the polarization operator is simpler. In this case we examine the response of the system 
to a total potential, and do not have to consider the additional gravitational forces arising 
from perturbations to the stellar orbits. Thus the Hamiltonian (|67D can be written as 

N 

Ho{Z) = lY.^^ + J2Mr,,ty, (79) 

i=l i 

here $j(r, t) = mj^j^i ^E'[r — rQj{t)] is the potential from all the stars other than i, moving 
along their unperturbed orbits. This Hamiltonian is separable, Hq{Z) = Y^iLi Hi{zi); the 
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responses of different stars are independent and hence we can work in 6- dimensional phase 
space instead of 6A^-dimensional phase space. There is a further simphfication in the hmit 
of large N, where the potentials $j can be replaced by the ensemble-averaged potential $0 
(the "mean-field" approximation, cf. §|^), and the analog to equation (^) is 

P(r, r', r) = e(r) J o!xdvFo(x, v)[5(x,(x, v) - r), 5(x - r')]. (80) 

Here Fo(x; v) is the one-particle df, defined so that Fo(X) v)(ix(iv is the (ensemble-averaged) 
mass in the phase-space volume rfxciv (in contrast to the A^-particle df, the integral of Fq 
over phase space is normalized to the total mass rather than to unity). The Poisson bracket 
is taken with respect to x, v; and x^(x, v) is the position at time r of the particle that was 
at (x, v) at time 0. 



4.1 Symmetries of the response operator 

For particular systems, the response and polarization operators may have symmetries in 
addition to those discussed in § pTT| . 

A common situation is that the equilibrium stellar system is invariant under time rever- 
sal (e.g. non- rotating galaxies). Under time reversal stellar positions remain the same, 
while velocities are reversed, (rj,Vj) — > (rj,— Vj). The density transforms as p{r,t) —>■ 
"^Sj'^(r ~ rj(— t)) = p(r, —t), while the Poisson bracket reverses sign. Consequently, under 
time reversal 

([p(r,t + r),p(r',t)]) ^ _([p(r, -t - r), p(r', -t)]) 

= -([p(r,t),p(r',t + r)]) 

= ([p(r',t + r),p(r,t)]). (81) 

Comparing to equation ([78|) , for systems invariant under time reversal we must have 

i?(r,r',r) = i?(r',r,r) (82) 

This result implies in turn that the Hermitian component of the response operator Rh{t^, i"', i^) 
is real and even in u, while the anti- Hermitian component Ra{t^, r', ^^) is imaginary and odd 
in uj. Similar relations can be proved for the polarization operator P(r, r', r), starting from 
equation (|80|) . 

Phase functions generally transform under time reversal as ■u[Z(t)] e^u[Z{—t)], where 
= +1 for mass or energy density, while = — 1 for momentum or angular momentum 
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density (e.g. Martin 1968). For systems invariant under time reversal, the general response 
function satisfies, 

i?™(r, r', r) = e^e^i?,„(r', r, r). (83) 
Likewise for stellar systems invarient under parity transformations (r^ —T^i), 

Ru.{r, r', r) = e^ef i?,,(-r, -r', r). (84) 

where is the signature of u{Z) under parity transformation. 



4.2 Goodman's stability criterion 



The results of the previous section can be used to re-derive an elegant instability test for a 
time-reversible stellar system (Goodman 1988). Consider the operator 

W(^) = P(^) - (85) 



where "^^^ = V^/(47rG) (c/. eq. ^5]). A collective mode 0s satisfies W(2;)(/)s = 0. Assume 
that the frequency z = is where s is real and positive. Then P{is) is real, and hence Hermi- 
tian if the stellar system is time-reversible. Similarly W is real and Hermitian. Therefore the 
eigenvalues A(s) of W(zs) are real. Moreover as s ^ oo, P(is) ^ so W — (27r)~-'^^^^, 
which is positive-definite for the inner product (^. Now let s decrease from infinity. If for 
any Sq > there is a function such that 

(0,W(zso)0) <0, (86) 

then W(zso) is no longer positive-definite. Thus some of its eigenvalues A(so) are negative. 
Therefore there is some Si > Sq for which one of the eigenvalues is zero. Hence there is an 
unstable mode with growth rate si. In other words a necessary condition for stability in 
time-reversible stellar systems is that 

27r(0,P(zso)0) < (0,*-V) (87) 



for all and for all sq > 0. This is a slight generalization of Goodman's result, which was 
derived only for systems with integrable potentials. 
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4.3 The polarization operator in action-angle variables 

The polarization operator ( ]80| ) can be evaluated explicitly if the potential $o(i") is regular, so 
that phase space can be described by action-angle coordinates (I, w). The orbits are given by 
I =constant, w = fit+wo where CI = dHo/dl, and mHQ is the Hamiltonian corresponding to 
$0- Jeans's theorem states that the one-particle df depends only on the actions, Fq = Fq{T). 
The canonical transformation from (x, v) to (I, w) conserves the volume element in phase 
space: (ixdv = (iwdl.^ 

Consider a single star of unit mass with action-angle coordinates (I, w) ; the correspond- 
ing spatial coordinate is x(I, w). Formally, the spatial density of the star can be written 

p(r|I,w) =5[r-x(I,w)]. (88) 

Since w is cyclic, we can expand this in a Fourier series, 

p(r|I, w) = ^pi(r|I)e^^-^ where pi(r|I) = — ^ f rfwe'^^'Xrll, w). (89) 

The Pi's are projection operators which give the Fourier components of any function of 
position g{r) that is expanded in a Fourier series of the form I]ifi'i(I) exp(il ■ w): 



-ilw 



^i(I) = J c^wc/[x(I,w)]e 

dwdr^(r)(5[r-x(I,w)]e-'''-^ 
dwdrg{r)p(r\I, w)e~*''^ 



{2ny 
1 

(2^ 
1 

(2^ 

drp,{r\l)g{r). (90) 



It is straightforward to prove the following useful identity: if h{I) is any function of the 
actions, 

(2^)^ E / dlh(l)p;{r\l)pi{r'\l) = 6{r - r') J d^fh[I{r, v)]. (91) 



1 

In particular, if h{I) is the one-particle df Fq{T), then 



(27r)3 Y: J rfIFo(I)K(r|I)p,(r'|I) = 6{r - r')po(r). (92) 



^We define the canonical momentum and the actions without the usual factor m, so these variables have 
dimensions (velocity) and (velocity) x (length) respectively. 
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The polarization operator (|80[) can now be written as 

P(r,r',r) = 0(r) j rfwdlFo(I) ^ [K(r|I)e"^i-(-+"^),p„,(r'|I)e*- 

l.m 

d 



= -{2nfeiT)^ J rfIFo(I)El- ^ [K(r|I)pi(r'|I)e~ 

= (2vr)3e(r)z /rflEl " ^^K(r|I)pi(r'|I)e-^'-"^ (93) 

where the last line follows through integration by parts. 
The Fourier transform of the polarization operator is 

where e is a small positive number. This can be split into Hermitian and anti-Hermitian 
components using the identity (|T^: 

P(r, r', u) = Puir, r', uj) + P^(r, r', u), (95) 

where (c/. eq. 

Ph{v,v\u;) ^ (27r)^^/dIp*(r|I)pi(r'|I)l-^p(j-^ 



dIK(r|I)p,(r'|I)l .-^6{l-n-uj). (96) 

Thus Pa is determined entirely by resonant stars satisfying 1 • O = o;. 

If the one-particle df depends only on the energy per unit mass E, Fq = Fo[E), then 
these expressions are simplified; for example, 

PAir,r\u;)=AnHujJ2 j dlpl{r\l)p,{v'\l)^5{\ ■ n - uj). (97) 

Equation (^) allows us to find an explicit expression for the rate of energy dissipation 
in a collective mode (eq. |63D , 

iy, = -47r^^^ j rfI|0,(I)|2l.^5(l.^^-^), (98) 

where the potential of the collective mode has been expanded as 0s(r) = X^i 0i(I) exp(il ■ w), 
and we have used (0) . This result is closely related to formulae originally given by Lynden- 
Bell & Kalnajs (1972) and can also be derived using second-order Lagrangian perturbation 
theory (Nelson & Tremaine 1995). Similar expressions to those in this subsection are also 
given by Tremaine & Weinberg (1984), Goodman (1988), and Palmer (1994). 



5 Fluctuations 



The density and potential of any stellar system fluctuate about their mean local values due to 
the finite number of stars. These fluctuations are described by the density-density correlation 
function and its Fourier transform. For an isothermal stellar system, it turns out that the 
correlation function is directly related to the response function. Thus, dissipation processes 
described by Ryi(cj) are determined by the fluctuation spectrum of the stellar system — the 
fluctuation-dissipation theorem. 



5.1 The correlation function 

We examine the density fluctuations Sp{r, t) in the stellar system, 

N 

Sp{r,t) = m^(5[r,(t) - r] - po(r), (99) 

i 

where po{r) — mN{S{ri — r)) is the mean density at r and as usual (•) denotes an ensemble 
average. The density fluctuations are characterized by the correlation function 

C{T,T',T) = {dpir,t + T)dpir',t)). (100) 

Since the equilibrium system is stationary, the correlation function is independent of time t. 
If we replace i by t — r we derive the symmetry relation 

C7(r,r',T)=C(r',r,-r). (101) 

The Fourier transform of the correlation function is called the dynamic form factor: 

S{r,r',u;) ^^J dre'^^C(v,v' ,r) ^^J dre'^^5p(r,t + r)5p(r',t)), (102) 

Alternatively, we can define S{r, r', a;) by expressing the correlation function in terms of the 
frequency transform of the density fluctuations, 

C(r, r', r) = y due'^^ j du'{5p{v, u)Sp{r', u;'))e"^^'"+'"'^*. (103) 

In order that the integral on the right side be independent of time t, we must have 

{5p{r, u})5p{r', -uj')) = S{v, r', uj)5{uj - J). (104) 
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The dynamic form factor satisfies the symmetry relations 

S*(r,r',cj) = S(r,r',-cj), S* {y,y\uj) = S{v' ,r,u)] (105) 
the first of these follows because the correlation function is real, and the second follows 



from ( [101| ) and implies that the dynamic form factor is Hermitian. Finally, the correlation 



function at zero time difference is given by integrating ^(r, r', uj) over all frequencies, 

C(r, r', 0) = y" duSir, r', uj) = (5p(r, t)6p{r', t)); (106) 

C(r, r', 0) is sometimes referred to as the static form factor. 

If the DF of the stellar system is invariant under time-reversal (c/. § [4 .21 ), and the Hamil- 
tonian has the form (p7|), the correlation function satisfies the principle of microscopic re- 
versibility, 

C(r, r', r) = C(r', r, r) or C(r, r', r) = C(r, r', -r). (107) 

This implies in turn that the dynamic form factor S'(r, r', u) is real, symmetric in r and r', 
and an even function of u. 

A simple proof of ( |107| ) relies on notation and results from We consider a generalized 
correlation function associated with any two phase functions Ms(Z), fs(Z), labeled by the 
parameter s. This may be written 

C™(s,s',r) = {5u,{t + T)5v,it)) = I dZU{Z)u,{Zr)v,,{Z) - {u,)^. (108) 



The density- density correlation function may be written in the form 

C(r,r',r) = j rfZ/o(Z)^;r(Z.)t;r'(Z) - Po(r)po(r'), (109) 

where fr(Z) = p{r,t) is defined by equation ([77|), and Po{r) = (t>r(Z)) is the ensemble- 
averaged density at r. Now change the dummy variable from Z to Z, where Z is obtained 
from Z by reversing all of the velocities. Since /o is time- reversible, we can replace /o(Z) by 
/o(Z). We can replace the volume element dZ by dZ because the Jacobian \d{Z)/d{Z)\ = 1. 
Moreover, fr'(Z) = fr'(Z) because depends only on coordinates, not momenta. Finally, 
(Z)t- = Z_T- because reversing the velocities yields the time-reversed orbit in the Hamiltonian 
(|67|); thus fr[(Z)^] = Vr{Z^r)- Equation ( p.09| ) becomes 

C{r, r', r) = J dZfo{Z)v,{Z.^W{Z) - po(r)po(r') = C(r, r', -r), (110) 



and equation ( |107| ) then follows from ( |101| ). 
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For some purposes it is useful to define a modified correlation function that removes 
long-term correlations: 

^(r, r', r) = C(r, r', r) - C(r, r', oo). (Ill) 

In contrast to most systems examined in statistical mechanics, the correlation function does 
not vanish at large times in most stellar systems. This is simple to understand in the 
context of spherical systems. The time-averaged density of each star is an annulus with fixed 
orientation and fixed inner and outer radii. The time-averaged total density is composed of 
the sum of the densities from N such annuli, and hence contains permanent irregularities 
that are not present in the ensemble-averaged density. We expect a non-zero correlation 
function at large times unless the system is ergodic. 

Correlations between particles are present for two conceptually distinct reasons: random 
fluctuations in the number of particles in a given small volume, which are present even if 
the particles move on their unperturbed trajectories; and gravitational interactions between 
particles, which deflect their mutual orbits. It is sometimes useful to isolate the correlation 
function arising only from random fluctuations in particle number, which we denote 

C(°)(r,r',r) = (5p(r, t + r)5p(r', t)) (112) 

As in equation (|80D , in the limit of large we can replace the exact potential $i by the 
ensemble- averaged potential $0) and write the correlation function as 

C(°)(r,r',r) =m J c^xdvFo(x, v)<5(x,(x, v) - r)<5(x - r') - lpo(r)po(r'). (113) 

The factor proportional to 1/N derives from a near-cancellation of terms analogous to the 
derivation of equation ( p.l8| ) below. 

Because the particle orbits are independent in this approximation, 6p{r, t) is a Gaussian 
random field, and hence its properties are completely described by the correlation function 
CW(r,r',r). 



5.2 The correlation function in action-angle variables 



We can find an explicit expression for the density-density correlation function C^°''(r, r', r) 
if the ensemble-averaged potential $o(r) is regular, using the same approximations and 
notation as in §4.3. 
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Assuming that the stars move on their unperturbed orbits in the potential $0^ we may 
write 



N N 



p(r,t) =m^p[r|I,,w,(t)] =m^^Pi(r|I,)e*'-("^*+^«). (114) 

1=1 i=l 1 

Since the initial phases Wjo are uniformly distributed, only the terms with 1 = survive 
in the ensemble average; thus the ensemble-averaged density is 

Po(r) = (p(r, t)) = {2Txf j dlF,{l)po{r\T) . (115) 

We may now write 

N 

5p{T,t + T)5p{r',t) = m^Y. EK(r|I*)pi'(r'|I,)e"''''"'^*^"^^"'°'^''''^"'*^"'"^ 

i,j=l 1,1' 

-p(r, t)po(r') - Po(r)p(r', t) + po(r)po(r'). (116) 

Now take an ensemble average; since the phases Wjo are randomly distributed, the average 
over exp[i(l' ■ wjq — 1 • Wjo)] will vanish unless either (i) i ^ j and 1 = 1' = 0, or (ii) i = j 
and 1 = 1'. In case (i) the sum over i ^ j yields A^(A^ — 1) terms, each of which is equal to 
Po(r)po(r')/iV2. Thus 

{Sp{r,t + T)Sp{T',t)) = ^^^K(r|I,)j9,(r'|I,) exp(-zl ■ ^r)^ - ^po(r)po(r'). (117) 

We can replace the sum over stars by the integral over the one-particle df to get 

C(o)(r,r',r) = {5pir,t + r)5pir' ,t)) 

= {2nfmJ2 J c?IFo(I)K(r|I)pi(r'|I)exp(-a-l^r)-lpo(r)po(r');(118) 

the superscript "0" is a reminder that we have neglected the self-gravity of the density 
fluctuations by using the unperturbed orbits of the stars. This result can also be derived 
from equations (|5), (|^) and (|113|) . 

The dynamic form factor ( |102| ) becomes 
5(°)(r, r', uo) = {2nfm J2 J rfIFo(I)K(r|I)pi(r'|I)5(l ■ Q - u) - lpo(r)po(r')5(cu). (119) 

Thus the dynamic form factor is determined entirely by resonant stars. 
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It is straightforward to calculate the static form factor ( |106|) , starting from equation 




(ixc?vFo(x, v)5(x — r)5(x — r') 



1 



Po(r)Po(r') 



N 



= m J (ixpo(x)(5(x — r)(5(x — r') — 
= mpo(r)5(r - r') - ^po(r)Po(r')- 



1 



Po(r)Po(r') 



N 



(120) 



This result also follows immediately from equations (0) and ( |118| ). 



The correlation function at large times is easily derived from equation (|118|) : 



CW(r,r',oo) = i2nfm J rfIFo(I)p* (r|I)j9o(r'|I) - -po(r)po(r'); 



(121) 



additional terms would appear if the potential had global resonances such that 1 ■ = for 
non-zero I (c/. Ranch & Tremaine 1996). Thus the modified correlation function defined in 



A system with Hamiltonian mHo^Z) is in thermal equilibrium if the A^-particle df has the 
form /o(Z) oc exp[— /3iJo(Z)] (in contrast the term "equilibrium" simply denotes that the 
one-particle df is a solution of the time-independent collisionless Boltzmann equation). For 
the usual gravitational force, the interaction potential m^\Ef(r) = — Gm^/|r|; in this case 
/o(Z) diverges exponentially as |rj — rj| — > and hence a thermal equilibrium state does not 
exist. Nevertheless, the useful concept of a thermal equilibrium can be retained for stellar 
systems, by modifying the interaction potential in one of two ways: 

• We can eliminate the interaction potential, setting \E'(r) = 0, so that the gravitational 
potential is determined entirely by the external potential $cxt(i'); we then augment 
$cxt(r) by adding the mean potential field of the unperturbed system. This "mean- 
field" approximation is equivalent to including the self-gravity of the equilibrium stellar 




(122) 



6 Thermal equilibrium and the fluctuation-dissipation theorem 
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system but neglecting the self-gravity of fluctuations, which is the same as approxima- 
tion (4) of §^ In this case the thermal equilibrium df is the product of one-particle 
isothermal dfs, 

/o(Z) (X Ul, exp[-/?i7(z,)], where H{z) = iv^ + $ext(r); (123) 

in other words the two-particle correlation function vanishes. This is a plausible model 
for galaxies, since (i) the one-particle df in galaxies is often approximately isothermal, 
and (ii) the relaxation times in galaxies are generally much larger than their age, so 
significant correlations between stars have not had time to develop. However, the mean- 
field approximation neglects the djTiamical effects of the self-gravity of the perturbed 
density; thus, for example, there is no distinction between the polarization and response 
operators and no collective modes. 

A more realistic approach is to soften the interaction potential to m^\l'(r) = —Gm^/ (r^+ 
^2^1/2 rjj_^g softening length h should be much smaller than the size of the system, so 
the large-scale equilibrium structure is unaffected; in fact it should also be much less 
than the typical interstellar separation so that softening does not affect most encoun- 
ters between stars. On the other hand the softening length should be large enough 
to suppress the strong correlations between stars that are present if the gravitational 
force is unsoftened. To see what is required, consider an infinite homogeneous stellar 
system, for which the equilibrium two-particle df is 

p(^)(zi,Z2) OC exp + + m^(ri - rs)]} 



P 



Zi)p('Hz2)[l + ^?(ri-r2)], (124) 



where 

p(^)(z) OC exp(-i/3t;2), ^((Ar) = exp [-/5m^(Ar)] (125) 

are the one-particle df and the two-particle correlation function. We would like the two- 
particle (and higher order) correlation functions to be small, so that the A^-particle df 
is the product of one-particle dfs as in equation (|123|) . This requires that max|/3m\E'| 



is small, which in turn requires b ^ Bq = Gm/a"^, where o"^ = 1/(3 is the mean- 
square velocity in one dimension. Fortunately, it is easy to satisfy these conditions in 
most stellar systems; thus, in the solar neighborhood we might (for example) choose 
b = 0.001 pc, which is much less than the typical interstellar separation of 1 pc but 
much greater than b^ ~ 10~^ pc. 



6.1 The fluctuation-dissipation theorem 
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We examine the relation between the response operator (|7^) and the correlation function 
( |108| ) when the stellar system is in thermal equilibrium. If /o(Z) oc exp[— /?iJo(Z)], 

= -/3/o(Z)^t;,,(Z_,); (126) 

in the second line, the replacement of the Poisson bracket with a time derivative follows from 
Hamilton's equations (|69D, and the argument of /o can be changed from Z_t- to Z because Hq 
and therefore /o is invariant along a trajectory. The response operator may now be written 

m (IT J m (IT J 

(127) 

in the second equality we have relabeled the variables Z and Z_^ as Z^ and Z respectively, 
and then changed the dummy variable from Z^ to Z; the Jacobian of the transformation is 
unity by Liouville's theorem. Finally, using equation ( |108| ), 

R^s, s', r) = ^^|-a.(s, s', r). (128) 
m OT 

This is the fluctuation-dissipation theorem (Callen & Welton 1951), which relates the re- 
sponse operator to the correlation function for any system in thermal equilibrium. 

If we set Us = Vs = p(s, t), where p(s, t) is defined by equation (ffTD, then the generalized 



correlation function Cuv and response operator Ruv of equations ( |108| ) and (^Tj) reduce to 
their analogs C and R defined earlier in equations ( |100| ) and (|^). Thus the fluctuation- 
dissipation theorem implies that 

i?(r,r',r) = ^^^C(r,r',r), (129) 
m OT 

or 

i?„(r, r', r) = :^-fc{r, r', r) = ;^|-C'(r, r', r), (130) 
2m OT 2m ot 

where Rq is defined in equation (p!3| ) and C is defined in equation ( |111| ). The formal analog 
of this result in Fourier space is 

RAir,r',u) = --^S{v,v',coy, (131) 

the equation is equally valid if ujS{r, r', u) is replaced by u!S{r, r', u) since the difference 
between the two functions is proportional to 006(00) which is zero. 



- 28 - 



If we suppress the interparticle gravitational interaction (the mean- field approximation), 
then the response operator R{r, r', r) is replaced by the polarization operator P(r, r', r), and 
the correlation function C(r, r', r) is replaced by C^^\r, r', r); thus 

P(r,r',r) = ^^C(°)(r,r',r). (132) 
m OT 

It is straightforward to verify this result when the ensemble-averaged potential is integrable, 
using the exphcit expressions for P(r,r',r) and C'^°^(r, r', r) in equations (p3D and ( |118| ). 
The analogous equation in Fourier space is 

PA(r,r',^) = -^5(°)(r,r',u;); (133) 

once again, the equation is equally valid if S^^'^ is replaced by S^^\ 



7 Correlations in nonisot hernial systems and the dressed-particle 

approximation 

Gilbert (1968) has described coUisional relaxation in stellar systems that are not necessarily 
in thermal equilibrium, by expanding the exact Liouville equation in powers of A^~^. His 
treatment is based on similar results derived by Rosenbluth et al. (1957) and Rostoker (1961) 
for a Coulomb plasma, and accounts fully for the self-gravity of the medium. 

These calculations show that fiuctuations in the system can be properly accounted for — 
to 0(A^^^) — by considering an individual test star and the polarization cloud it induces in 
the background medium together, as a "dressed particle" ; for a test star with trajectory r*(t) 
the density of the bare particle is po{'r,t) = m6[r — r^{t)], and the density of the dressed 
particle is given by 

Prf(r, uj) = [I + 27rR(cu)*]po(r, oo) = A(cu)po(r, oo) (134) 

(c/. eq. The test star itself is drawn from the one-particle df and moves on its un- 

perturbed orbit; thus in this approximation the only fluctuations in the density of dressed 
particles are statistical fluctuations due to the flnite number of particles. The papers de- 
scribed above effectively show that induced correlations between dressed particles are 0(A^~^) 
or higher. 
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In the dressed-particle approximation, the dynamic form factor S can be computed from 
S^^-*, simply by replacing each particle by a dressed particle: starting from equation ( |104| ) 

S{uj)6{uj - uj') = {6pd{r,uj)6pd{r',-uj')) 

= AHS(°H^)A^(cu)5(cu-cu'). (135) 

The advantage of this replacement is that S'-''-* can be computed directly from the unperturbed 
orbits of the particles (e.g. eq. |119| , if the ensemble-averaged potential is regular). 

These results are consistent with the two forms of the fluctuation-dissipation theorem 
involving the response and polarization operators (eqs. |131| and |133|) , in the following sense. 
The latter equation can be written as cuS^^^ = {2mi/P)PA', substituting this into equation 
( |135|) gives 

= ^AP^At. (136) 

Using equation (]26|) we can show that APA^ = R+27rR*R^ so that APaA"^ = (APA"^)^ = 
Ha] hence equation (|136|) simplifies to 



^ 2mi 

ujS = ^Ra, (137) 



which is the same as equation ( |131|) . 



Ichimaru (1965) argues that the result (|135| ) should hold for any stationary plasma, 
whether or not in thermal equilibrium, so long as the relaxation time is much longer than 
the correlation time; in the context of stellar systems, the analogous constraint is that the 
relaxation time should be much longer than the crossing time, which is equivalent to the 
modest requirement that the number of stars ^ 1 and entirely compatible with the results 
of Gilbert (1968). 



7.1 Fluctuations and collective modes 



One consequence of equation ( |135D is that the density fluctuations described by S{u) become 
very large for frequencies near singularities of A(z) or Ii.{z), i.e. near the frequencies of 
collective modes. To see this more explicitly, we can expand H{z) in a Laurent series near a 
collective mode with complex eigenfrequency zo = ujo + if]- The resonant part is then 

R,es(^) = (138) 

Z — Zq 
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where R_i is the residual of H a,t zq. Consequently, the dynamic form factor at real fre- 
quencies near a weakly damped collective mode is 

{uj-uJor + T]^ \t]\ 
where A_i = 27rR_i^', and we have used equation (|T 



Note that the level of fluctuations becomes very large as rj —>■ — that is, as the mode 
approaches neutral stability. This phenomenon is analogous to opalescence in the vicinity of 
a critical point. Large fluctuations have been observed in numerical simulations of systems 
near marginal stability (Ivanov 1992, Weinberg 1993). 



8 Applications of the fluctuation-dissipation theorem 

A useful preliminary calculation is the response of a stellar system in thermal equilibrium 
to a static or slowly growing external potential $e(r). Using equations (|l|), ( |129| ) and (|111|) 
we find that the density perturbation induced by $e is 

//•oo 
dr' rfri?(r,r',r)<l>e(r') 

m J Jo OT 

■ J dr%e{r')C{r,r',0). (140) 



A 

m 



If, in addition, we ignore induced fluctuations and approximate (5(r, r', r) by C'*^°^(r, r', r), 
and the ensemble-averaged potential is integrable, we find using equations ( |122| ) and ( PU[ ) 
that 

pf(r) = -(27r)3/5^/rfIFo(I)K(r|I)$e,i(I) 

= -/? 1 rfvFo(r, v) [<l>e(r) - ; (141) 

here ($e) = *^'e,o[I(i", v)] is the orbit-averaged potential experienced by the particle passing 
through the phase-space point (r,v). The last line is a special case of the general result 
(e.g. Lynden-Bell 1969) that the linear response of a stellar system whose df depends only 
on energy E to a, slowly varying potential is 

/=^(4-(4)). (142) 
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8.1 Dynamical friction on an orbiting body 



Consider a body whose center of mass travels on an orbit r^,(t) through a stellar system, and 
whose gravitational potential is $*(r — r*) (normally is spherically symmetric but this is 
not necessary for the derivation; we do, however, neglect changes in the orientation of the 
body). We wish to compute the mean force Fi(t) exerted on the body by the response it 
induces in the stellar system. 

The density in the wake, ps{r,t), is given by equation (|l]), and Fi is equal and opposite 
to the force exerted on the wake by the body, that is. 



Fi(t) = / drV$*[r- r,(t)]p,(r,t) = / drc?r'c?t'V$,[r - r,(t)]/?(r, r', t - t')'^'*[r' - r,(t')] 



For completeness we point out that there is another component of the mean force on 
the body, which is not caused by the response it induces in the system. The fluctuating force 
on the body is given by 



Although {6p{r,t)) = 0, the mean fluctuating force does not vanish because the orbit r^:(t) 
is correlated with the fluctuations 6p{r,t); thus there is an additional contribution to the 
mean force F2 = {5F). 

To compute F2, we flrst compute the perturbation to the orbit caused by the fluctuating 
field, r*(t) = ro(t) + Sr^(t) where 




(143) 




(144) 





r) (145) 



plus higher-order terms. 



Expanding equation 



( [L44| ) and integrating by parts gives 




(146) 



where 



{*r.(()-Vp(r,()> 




Jo 



00 



drrV'^^r' - r^t - r)] ■ V(5p(r, t)6p{r', t - r)) 




(147) 
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Thus, 

F2(t) = -^ /drcir'V$,[r-r,(t)] T drrVC(r, r', r) • V'<l>4r' - r,(t - r)]. (148) 
M J Jo 

Note that F2 can be formally divergent if C(r, r', 00) is non-zero. The reason is that in this 
case the time-averaged density is not the same as the ensemble-averaged density, so that the 
"fluctuating" force can contain a constant component that leads to a secular change in the 
orbit. 

The two components of the mean force, Fi and F2, have different dependences on the 
mass of the orbiting body; in particular, if this mass is much greater than the mass of the 
stars in the system, M 3> m, then the dominant force is Fi, which is greater than F2 by 
0(M/m). In the opposite limit, M — ^ and oc M, the mean acceleration ai = Fi/M 
goes to zero as M while the acceleration a2 remains finite. 

If the stellar system is in thermal equilibrium, we can use the fluctuation-dissipation 
theorem in the form (|129|) to write equation (|143|) as 

Fi(t) = ^ / rfrt/r'V$4r-r,(t)] T rfr^^^^%^<l>,[r' - r,(t - r)]. (149) 
m J Jo or 

Now replace C(r, r', r) by C(r, r', r) = C(r, r', r) — C{r, r', 00) and integrate by parts with 
respect to r: 

Fi(t) = / rfrrfr'V$4r - r,(t)]C(r, r', 0)$4r' - r,{t)] 

I drdr'V^,[v-Y,{t)] r dTC{r,r',T)v,{t - t) ■ V^,[r' - r,{t - t)] 

TYXi tj '-'0 

= F, + Frf, (150) 

where v^, = dr^^/dt. Using equations ( |14CI| ) and ( |144| ) it is easy to show that the first term, 
Fs, is simply the force on the body due to the static response of the stellar system; in other 
words this is the force that would result if the body were fixed in its present position. The 
second term, F^^, vanishes if v^, = and corresponds to dynamical friction. The dynamical 
friction force can be rewritten as 

FdAi) = -- r dTv„{t-r){6F,{t)6F^{t-r)), (151) 
m Jo 

where the fluctuating forces 5Fi{t) on the body are computed using the modified correlation 
function C(r, r', r). The result ( [151| ) holds for any stellar system in thermal equilibrium. 

In most calculations of dynamical friction the self-gravity of the wake is neglected. 
In this approximation the derivation of equation ( |151| ) remains the same, except that the 
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fluctuation-dissipation theorem is used in the form ( p.32| ) and the fluctuating quantities 5p 
and are computed using the correlation function C^^\ which neglects interactions between 
stars. 



8.2 Infinite homogeneous medium 



It is instructive to work out the mean force and fluctuating force on a point mass M traveling 
through an infinite homogeneous system of stars of mass m. We neglect the self-gravity of the 
equilibrium system, so that all objects travel on straight-line orbits at constant velocity; we 
also neglect the self-gravity of the response, since otherwise the equilibrium system is Jeans 
unstable. Thus there is no distinction between the polarization and response functions P 
and i?, or between the correlation functions C^^^ and C . 

Following equation (173), the response function may be written 



i?(r,r',r) = ^([p(r, t + r), p(r', t)]) 



m 



W E([^(r - r.o - v,or), 6{v' - r,o)]), (152) 



where the orbit of particle i is rj(t + r) = r-jo + VjoT. The Poisson bracket is taken with 
respect to the phase-space variables (fjo, Vjo). Thus 

[6{r - r,o - v,or), 6{r' - r,o)] = -6,,rV6[r - r,{t + r)] ■ V'6[r' - r,{t)]. (153) 

The response function is then 

R(r, r', r) = -me(r)r ^ {V6[r - nit + r)] ■ V'6[r - rj{t)]) = -^^tV ■ V'C(r, r', r). 



00 

m 

(154) 



Integrating equation ( |14^ ) by parts with respect to r', substituting equation (|154|) , and 
comparing to (|143| ) shows that in an infinite homogeneous medium there is a simple relation 
between the components of the mean force, 

F.^(^)F. (155) 

Thus 

F = Fi + F2= (^1 + ^) y'rfrrfr'^"drV$4r-r,(t)]i?(r,r',r)<I>,[r'-r,(t-r)]. (156) 
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To proceed further, we shall assume that the stellar system has a uniform number 
density n and a Maxwellian velocity distribution with one-dimensional dispersion a — which 
therefore is isothermal, with /? = 1/cr^. In this case the force correlation function can be 
shown to be (c/. Cohen 1975) 



(Fi(r,t + r)Fj(r',t)) =Txn{GmMfx 

, . 2 



-,-2) eriiu) - ^e-«^ 









a 





-e 



(157) 



where a = r' — r, m = a/(2^/^cr|r|), and erf(a;) = 271^^^"^ Jq exi){—y'^)dy is the usual error 
function. The trace of ( |157| ) is simpler, 

(F(r, t + r). F(r', t)) = '"''^'"^''% rf(.). (158) 

d 

an expression given by Cohen (1975). 

For a homogeneous system the static induced force F^ (eq. p. 501 ) vanishes, and we have 

F, = (l + ^) F,, = - (l + ^) £^ £ dr{5F,{r,t + r)5F,(r - v.r,t)); (159) 

we have changed the limits of integration to ri > and r2 < oo to avoid divergences. The 
integral is easily evaluated to yield 

— A-KG^MmiM + m)n\nA ^ ^, , 

F = -V, ^— [erf (m) - u erf'(M)] , (160) 

^* 

where u = f*/(2^/^cr) and A = T2/T1. This is the standard formula for dynamical friction 
in an infinite medium with Maxwellian velocity dispersion (Chandrasekhar 1943, Binney & 
Tremaine 1987), except that the Coulomb logarithm In A is determined by the limits on 
the time lag in the correlation function rather than by limits on the impact parameter. 
Note also that the frictional force induced by the body — which is what we defined to be 
dynamical friction in the previous subsection — is actually Fi, not F, which is smaller by a 
factor M/(M + m). 

The effects of relaxation in an infinite homogeneous medium are often expressed in 
terms of diffusion coefficients, (Afj) and (AviAvj), representing the mean and mean-square 
changes in velocity per unit time At. These changes are related to the forces on a point 
mass M by 

/■At i-At i-At 

M{Avj) = / Fj{t)dt, M^Av^Avj) = dt dt' {5Fi{t)5FM')). (161) 

JQ Jo Jo 
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If we assume that the correlation function is neghgible for time lag > Tmax? where Tj^^x ^ 
At — this is the Markov approximation of §1 — then the second of these equations may be 
written 



M^{AviAvi)= dt (5FAt)5FAt - t)) = 2 dt (5FAt)6FAt - t)), (162) 

"'0 ^-rn,ax Jo Jo 

where the second equation follows from (|107|) . Equation (|159|) then yields 



{Av,) = -v„{Av,Av,), 163 



a result derived already by Chandrasekhar (1943). 



9 Discussion 



We have described stellar dynamics using the language of statistical physics: linear response 
operators, correlation functions, and the fluctuation-dissipation theorem. It is fair to ask 
what we have gained from this rather formal approach, especially since stellar systems do 
not satisfy many of the simplifying assumptions commonly used in many-body theory, such 
as homogeneity, local forces, and thermodynamic equilibrium. We believe that there are two 
main reasons why pursuing this approach is worthwhile. First, these techniques have proven 
to be extremely powerful in other branches of physics, so it is important to understand to 
what extent they can be applied to gravitating A^-body systems. Second, it is useful to know 
which features in stellar systems result from general properties such as causality and time- 
reversal symmetry, and which depend on specific approximations used to make the dynamics 
tractable. 

An occasional controversy, for example, is the appropriate maximum impact parameter 
that should be used in the Coulomb logarithm that enters the diffusion tensor computed 
using the local and Markov approximations of §[1|. Chandrasekhar (1942) and Kandrup 
(1981) advocate terminating the integration over impact parameters at the typical interstellar 
separation, whereas Cohen et al. (1950), followed by most modern authors, argue that 
the integration should include all impact parameters up to the characteristic size of the 
stellar system (or the Debye length in a plasma). In an iV-star system with N ^ 1, the 
polarization and response operators depend only on the one-body df, which is unchanged if 
the number of stars is changed so long as the system mass Nm is conserved. The response 
operator determines dynamical friction, and the diffusion tensor is related to dynamical 
friction through ( |151| ). Thus the diffusion tensor is unchanged — except for small terms 
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whose fractional contribution is 0(m/M) — if the number of stars is changed, so long as Nm is 
conserved. However, this change affects the interstellar separation; hence the diffusion tensor 
cannot depend on the interstellar separation and the effective maximum impact parameter 
must be of order the system size. 



A central result of this paper is equation (|151|) , which relates the dynamical friction force 



to the force-force correlation function in any stellar system in thermal equilibrium. This quite 
general result illuminates the relation between stochastic and dissipative gravitational forces 
in stellar systems. For example. Ranch & Tremaine (1996) have argued that the rate of 
angular momentum relaxation is strongly enhanced in nearly Keplerian star clusters, such as 
those found around massive black holes ( "resonant relaxation" ) . Equation ( |151| ) immediately 



implies that dynamical friction is similarly enhanced in such clusters, an effect analyzed by 
Ranch and Tremaine and termed "resonant friction" . 

We have not discussed more general issues related to the long-term relaxation of the 
one-particle df, as described by the full collisional Boltzmann equation or its generaliza- 
tions. For isothermal systems, the fluctuation-dissipation theorem already suggests that the 
Boltzmann collision integral describing this relaxation must be related to the power spec- 
trum of potential fluctuations in the background medium. Indeed, in plasma physics, the 
corresponding Balescu-Lenard collision term can be expressed in terms of the fluctuation 
spectrum determined by the collisionless plasma (e.g. Lifshitz & Pitaevskii 1981). Weinberg 
(1993) has derived this collision integral for a model periodic system, and argues that fluctu- 
ations associated with nearly unstable collective modes can strongly enhance the relaxation 
process. 

The concepts in this paper can equally be applied to study artificial fluctuations and 
dissipation in numerical methods that approximate stellar systems, such as self-consistent 
field codes (Hernquist & Ostriker 1992), which approximate the gravitational field by a 
truncated multipole expansion. 

The fluctuation-dissipation theorem in statistical physics also relates transport coeffi- 
cients (e.g. electrical and thermal conductivity, diffusion coefficients) to frequency moments 
of the fluctuation spectrum. In many cases, these transport coefficients can be calculated 
without full knowledge of the correlation function. Moreover, conjugate transport coefficients 
are related by the Onsager relations. Similar relationships may exist for stellar dynamical 
systems, although their usefulness remains unclear, especially for non-isothermal systems. 
To study the dynamics of a stellar system described by a general one-particle df, it may be 
possible, for example, to consider phenomenological response operators that are consistent 
with the lower-order moments of the exact response function. These models may offer an- 
alytic insight into the linearized dynamics of the stellar system without computing the full 



perturbed distribution function. 



10 Summary 



The main goal of this paper has been to assemble and discuss the general properties of linear 
response, dissipation and fluctuations in stationary stellar systems. Many of the results are 
not new, having already been derived in other arenas of statistical mechanics — although 
usually in the context of spatially homogeneous systems with short-range forces. Here we 
summarize our main results. 

We have expressed the dynamical response of a stellar system in terms of a linear re- 
sponse operator R{r, r', r), which determines the density perturbation induced by an external 
potential, 

p,(r,t) = J dr'dt'R{y,v\t-t')^,{v\t'). (0) 

The response operator is analogous to the conductivity tensor used in plasma physics; in 
statistical physics it is sometimes called the generalized susceptibility (Landau & Lifshitz 
1980). Many of its analytic properties in the frequency domain follow from the requirement 
that the response must be causal. For example, its Hermitian and anti-Hermitian parts 
Rj^(a;) and RA(t^) are related through the Kramers-Kronig relations (p!9|). 

The anti-Hermitian part of the response operator is associated with dissipation and 
dynamical friction; it determines, for example, the irreversible work done on the system 
by an external time-dependent potential (§ p.3|) . If the perturbation is periodic, (^eij^.t) = 
Re[0e(r)e~"^°*] the rate of irreversible work done is 

W = iriujo J dr(ir'0*(r)i?^(r,r',cJo)0e(r') = vr2CJo(0e, R-A^e)- 

For a Hamiltonian system, an exact formal expression for the response operator is given 
by the ensemble average of a Poisson bracket in 6A^-dimensional phase space, 

i?(r,r',r) = ^([p(r,t + r),p(r',t)]), 
m 

where p{r,t) = mX]i^[r — is the exact density distribution of the system. 



Equation (ffSD generally cannot be explicitly evaluated for realistic stellar systems (al- 
though short-time expansions and sum rules can be calculated). However, a closely related 
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polarization operator P{r, r', r) is easier to evaluate. The polarization operator relates the 
induced density perturbation to the total (external plus induced) potential, 

If the ensemble-averaged potential of the stellar system is integrable — so that individual 
stellar orbits have well-defined actions I — the polarization operator is given by 

,p;(r|I)p,(r'|I) , dF„ 



P(r,r',.) = (2.)=E/<'lf«^l- 

,j I ■ il — te — uj 



dl 



where -fo(I) is the one-particle df, and the functions pi(r|I) are projection operators onto 
action space (eq. 



The operators R and P are related by a nonlinear equation (^61). The operator P 
determines the dispersion relation for collective modes of the stellar system (^61). The rate 
of energy loss of a mode to Landau damping is determined by its anti-Hermitian part, 

Wm = vriicj(0„ Pa(Ps) = 2//^™, (ID 

where and t] are the energy and damping rate of the mode. Dissipation occurs through 
resonant interaction with stellar orbits commensurate with the mode frequency. 

Density fluctuations in a stellar system are characterized by the correlation function 

C(r,r',r) = (5p(r, t + r) V(r', t)). 

For a stellar system in thermal equilibrium, described by an A^-particle distribution function 
/o(Z) oc exp[—/3Ho{Z)], the response function is directly related to the correlation function 
by the fluctuation-dissipation theorem, which states that 

i?(r,r',r) = ^^|-C(r,r',r). O 
m OT 

The equivalent expression in the frequency domain reads 

i?^(r,r',cu) = -^5(r,r',^), O 



where S is the Fourier transform of the correlation function (eq. |102| ). Thus dissipational 
processes in the medium, determined by R^(lo'), are directly related to the power spectrum 
of density fluctuations. 
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